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ABSTRACT 
Structural  design  problems  can  be  considered  to  be  optimization 
problems  because  a  design  is  sought  which  is  optimal  by  some  criterion 
subject  to  limitations  on    size,  behavior,  or  other  aspects  of  the  struc- 
ture.   Under  certain  conditions,  such  problems  may  be  solved  by  con- 
ventional mathematical  programming  techniques.    The  minimum  weight 
design  of  a  statically  indeterminate  three-bar  truss  is  used  to  illustrate 
the  application  of  the  "sequential  unconstrained  minimization  technique" 
of  Fiacco  and  McCormick  to  the  optimal  design  of  trussed  structures.    A 
suggested  computational  procedure  and  computer  program  are  included. 
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I.     INTRODUCTION 

The  desired  conclusion  of  the  structural  design  process  is  a  structur- 
al configuration  which  performs  its  intended  function  efficiently.     The 
method  by  which  the  final  design  is  determined  usually  includes  the 
establishment  of  an  initial  configuration  from  which  the  final  design  is 
evolved  through  a  process  of  analysis  and  redesign. 

The  electronic  computer  has  permitted  an  acceleration  of  the  design 
process  by  enabling  the  rapid  analysis  of  the  mathematical  models  assoc- 
iated with  the  structures  under  consideration.     In  order  to  more  efficiently 
use  this  capability  L.  Schmit  [1]   has  suggested  a  rational  approach  to  the 
design  process  called  "systematic  structural  synthesis."    This  approach 
involves  the  systematic  evaluation  and  modification  of  a  large  number  of 
trial  designs  until  the  optimal  design  is  obtained. 

The  basic  problem  is  to  design  a  structure  which  is  optimal  by  some 

criterion  subject  to  a  set  of  requirements  which  specify  acceptable  limits 

on  the  behavior,  size,  or  other  aspects  of  the  structure.    Such  problems 

are  known  to  the  operations  analyst  as  mathematical  programming  problems, 

the  general  form  of  which  is  as  follows.     Find  a  vector  X  which 

minimizes    F(X) , 

subject  to    Hi(X)    -    0,  i  =  1 ,  2  ,  .  .  .  ,m; 

Gj(X)     £    0,  j  =  1,2 ,p. 

T 
Here,    X    =    (x1,x2,  ...,x  )      is  an  n-dimensional  column  vector; 

F(X)  is  the  objective  function;  and  the  relations  Ui  (X)  =  0  and  Gj(X)  ^0 

represent  the  constraints  on  the  problem  s 


If  F(X)  is  convex  in    X    and  the  constraint  region  is  convex,  i.e.  the 
set  of  all  solutions  to  H±  (X)  =  0 ,  i  =  1 , .  .  .  ,m  and  Gj  (X)  =  0 ,  j  =  1  ,  .  . .  ,  p 
forms  a  convex  set,  the  problem  is  called  a  convex  programming  problem. 
If  a  problem  can  be  classified  as  a  convex  programming  problem  then 
there  are  solution  procedures  which  are  guaranteed  to  find  the  global 
optimum. 

Several  techniques  which  have  been  devised  to  solve  this  special 
case  of  the  general  programming  problem  are  Rosen's  "gradient  projection 
method,"  Zoutendijk's  "method  of  feasible  directions,"    and  Kelly's 
"cutting-plane"  method.    These  procedures  are  summarized  by  Hadley  in 
Ref.  7.    Another  convex  programming  technique  called  the  "method  of 
alternate  steps"  was  suggested  by  Schmit  [l]   as  a  means  for  converging 
on  the  optimal  structural  design. 

A  recently  developed  convex  programming  technique  which  appears 
to  be  especially  promising  is  the  "sequential  unconstrained  minimization 
technique  for  convex  programming  with  equality  constraints"  of  Fiacco 
and  McCormick  [2]  .    By  this  method,  the  objective  function  and  con- 
straints are  dealt  with  simultaneously  by  combining  them  into  a  single 
function.    This  function  is  minimized  for  different  values  of  an  auxiliary 
parameter  thereby  generating  a  sequence  of  feasible  solutions  that  con- 
verge to  the  optimal  solution. 

If  the  convexity  restrictions  for  the  constraint  region  or  the  objective 
function  are  relaxed  a  convex  programming  solution  technique  may  still  be 
used  but  there  is  no  guarantee  that  a  solution  to  the  problem  is  other  than 
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a  local  minimum.    For  most  practical  problems,  it  is  no  less  valuable  to 
obtain  information  about  local  optima  in  the  absence  of  a  global  solution. 

In  general,  structural  design  problems  do  not  meet  the  convexity 
restrictions  for  convex  programming  problems.    One  such  problem  is  the 
determination  of  the  cross-sectional  areas  of  the  members  of  a  statically 
indeterminant,  coplanar,  three-bar  truss  such  that  the  total  weight  of  the 
truss  is  minimized  subject  to  constraints  on  stress  and  displacement. 
This  structure  is  illustrated  in  Figure  1.     Procedures  which  can  be  used 
to  solve  this  problem  can  also  be  applied  to  other  trussed  structures. 
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II.     OBJECTIVE  AND  SCOPE 

The  purpose  of  this  study  was  to  illustrate  the  application  of  the 
Fiacco-McCormick  technique  to  the  optimal  design  of  trussed  structures. 
Since  minimum  weight  has  usually  been  the  major  goal  in  the  development 
of  aircraft  structural  design  techniques,  it  was  used  as  the  criterion  for 
optimality . 

One  example  of  the  determination  of  the  minimum  weight  design  of  a 
three-bar  truss  is  included  with  results  and  the  computer  program.    Also 
included  is  a  brief  discussion  and  outline  of  the  solution  procedure  used 
to  solve  the  problem.    This  procedure  can  only  be  applied  to  problems  for 
which  constraints  can  be  formulated  in  completely  analytical  form. 
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HI.     THE  THREE-BAR  TRUSS  PROBLEM 
The  three-bar  truss  problem  of  Schmit  [l]is  illustrated  in  the  following 
sketch.  y 


Figure  1.    The  Three-Bar  Truss 


The  function  of  the  truss  is  to  transmit  the  load  Pj  from  the  point 

of  application  at  joint  A  to  the  support  B.    The  deflections  of  the  joint 

A  in  the  x  and  y  directions  due  to  the  load  P,  are  u   .  and  u    .  respectively. 

J  xj  yj 

For  a  coplanar  system  of  forces ,  one  of  the  conditions  for    equili- 
brium of  a  rigid  body  in  one  plane  is  that  the  algebraic  sum  of  the 
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components  of  all  forces  in  each  of  two  mutually  perpendicular  directions 
in  the  plane  of  the  forces  must  separately  vanish.    Therefore, 

£F       =     P-lj      COS^S-l   +   p2j    COS02  +   p3j      COS/33      +   Pj    COSftj    =   0 

and 

EF     =  pXj  sin/Sj  +  p2j  sin^s  +  p3J  sinjS3  -  Pj  sinc^  =  0 
where  ptJ    is  the  force  component  in  the  direction  of  the  ith  structural 
member  due  to  the  jth  load  Pj  . 

The  stress  in  the  ith  member  due  to  the  jth  force  is  defined  as 

°U  =  Pij  A 
where  Aj    is  the  cross-sectional  area  of  the  ith  member.    The  equili- 
brium equations  can  then  be  written  as 

A1o1J   cos/3-l  +  Ajgcrsj  cosjS2  +  A3  a3j  cos/S3  +  Pi  cos&j  =  0  (1) 

and 

Aloli  sin^  +  A2a2j  sinj52  +  A3o3j  sinj33  -  Pj  sinttj  =  0  (2) 

Assuming  small  displacements,  the  change  in  length  of  the  ith 

structural  member  due  to  the  jth  load  is 

6  i,  =  -  u    .  cosfi  -  u    .  sinpi   .  (3) 

*j  xj  yj 

A  relationship  between  stess,  displacement,  and  temperature  change  can 
be  established  if  is  noted  that  6^  may  also  be  expressed  as 

6'J  =  aT^T    +   ^i^T,  lt  .  (4) 

where  lt  is  the  length  of  the  ith  member,  a  i  is  the  mean  coefficient  of 
thermal     expansion  of  the  ith  member,    Et  is  the  modulus  of  elasticity  of 
the  ith  member,  and   ^rYi  is  the  temperature  change  imposed  externally  on 
the  ith  member. 
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If  the  temperature  is  assumed  to  be  constant,  equations  (3)  and  (4) 
can  be  combined  to  yield 

Pijli  . 

— — —    +  u    .  cosp.  +  u    .  sin/3,  =  0 
Ai  Ej  xj  »         uj         ^ 

and  since  ali  =  P^/Aj  and  lt  =  N/sin/Sj  ,  this  equation  becomes 

„       .    a        +  u    .  cosjS,  +  u    .  sinjS,  =  0.  (5) 

Ej  sinPi  xj  yj 

For  purposes  of  this  illustration,  it  is  assumed  that  all  geometric 

parameters  of  the  truss  are  given  except  for  the  cross-sectional  areas. 

It  is  further  assumed  that  in  addition  to  non-negativity  restrictions  for 

the  At  values,  the  behavioral  limits  for  the  structure  are  specified  by 

upper  and  lower  bounds  on  the  stresses  and  displacements.    Then,  for 

the  jth  load  condition, 

LU  *  °a  ^UtJ  ,  i  =  1,2,3, 

L   .  <  u        ^  U    ., 
xj         xj  xj 

and 

L    .  <  u       £  U   .. 

yj       yj        yj 

If   pt  denotes  the  specified  density  of  the  material  in  the  ith  struc- 
tural member  then  the  weight,   W,  of  the  truss  is 

w  =  Ai  la  Pa  +  A2  12  P2  +  A3  13  p3 
Since  li  =  INj/sin/^  ,  this  equation   can  be  rewritten  as 

w    =     Ai    N  Pi     +      A2    N   p2      +  A3   N   p3 
sin^a  sin£2  sin£3 

which  is  linear  in  A  and  is  therefore  convex. 
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The  mathematical  programming  problem  is  then  to  fine  Al  ,  A2  ,  A3  , 

the  cross -sectional  areas  of  the  structural  members  which  minimize 

N*  p,  N*  pz  N-Oo 

W       =      A,       ,     s  l   +    A2        ,    a      +  A3       ,    ■§ 

1    sin^  sinp2  sinp3 

subject  to  the  constraints 

Ai°ij  cos^  +  AgUaj  cosi?2  +  A3a3J  cos/>3  +  Pj  coso^  =  0 

Aloli  sinp\  +  A2or2j  sin/32  +  A3a3j  sin/33  -  P^  sinc^j   =  0 

N 

a,,     -r —  +u    .cosjSt  +  u    .sine,        =0 

1J     E1  smjSi  xj        ^        xj         x 

°2>     E8sinft,  +  ux.cos^2  +  uy.sinfe      =0 

as'  Easinfc         +  "xj003^  +  "yjsin&       =  0 
and 

Lutein    *   Ui,  ,  L   ,    s   u        £  U    .,    L        <    u        <s    U    .,  A,  ^0 
1J         1J  1J         xj  xj  xj         yj  yj  yj 

for  i  =  1,2,3  and  j  =  1 ,2  ,  .  .  .  ,k  where  k  denotes  the  number  of  distinct 
loading  conditions.    It  is  assumed  that  the  loading  conditions  are  not 
applied  simultaneously.     Note  that  inequalities  such  as  LtJ    ^  atJ    <,  ut. 
can  be  rewritten  as 

Git     -     tu     *    0, 

-ati   +    u4J   a  0. 
Now,  in  this  problem  statement,  the  variables  to  be  determined  are 
actually  A1  ,A2  ^3,0^  ,  1=1,2  ,3,  j  =  l ,  .  .  .  ,k  and  u    .  /  u    .,j  =  l,...,k. 
It  is  therefore  apparent  that  the  first    two  equality  constraints  are  non- 
linear because  of  the  cross  products  o^  At  and  hence  the  constraint 
region  is  non-convex. 
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IV.     SOLUTION  PROCEDURE 
Recall  that  the  general  form  of  the  mathematical  programming  problem 
is  to  determine  the  vector  X  which 
minimizes      F(X), 
subject  to      Hj  (X)  =  0 ,  i  =  1 ,2  ,  .  .  .  ,m; 

Gj  (X)  *0,    j  =  1,2, p, 

T  T 

where  X    =  (xj  ,x3  , .  .  .  ,xn)      is  an  n-dimensional  column  vector. 

The  basis  of  the  Fiacco  and  McCormick  technique  is  the  P-function 

which  is  defined  as  follows  for  this  general  form. 


m 


1  -     -1  ■  *a 


P(X,r)  =F(X)  +r       £      Q    (x)       +  r*     I  (H,  (X)  )■  . 

i=l  j=l 

The  general  solution  procedure  is  then  to  seek  minimum  values  of  P(X,r) 

for  r  =  t1  ,r2  ,  .  .  .  ,rk    where  rl  >r2  . .  .  rk  >  0  .    Then 

lim  Xk  =  X 

and 

lim  P(Xk   ,rk  )  =  F 
k-->  °° 

where  F    is  the  minimum  value  of  the  objective  function  of  the  mathema- 
tical programming  problem  at  X. 

The  conditions  for  the  determination  of  a  global  minimum  for  the  P 
function  are  generally  the  same  as  for  the  convex  programming  problem. 
The  objective  function,  F(X),  must  be  a  convex  function.    Recall  that  a 
linear  function  is  convex  although  not  strictly  convex.     Also  the 
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constraint  region  must  be  convex  (which  is  not  true  for  this  problem) 
and  F ,  G1 ,  .  .  .  / GB ,  Hl ,  .  .  .  ,  H     must  be  continuous  . 

Close  inspection  of  the  P  function  reveals  that  the  term 
m 


r      I 


Gi(X) 

i=l 

is  a  boundary  repulsion  term  which  keeps  the  minimum  feasible  solution 
of  the  P-function  in  the  interior  of  the  region  defined  by  the  inequality 
constraints.     Naturally,  the  use  of  such  a  boundary  repulsion  term 
requires  that  solutions  exist  in  the  interior  of  the  region  defined  by  the 
inequality  constraints. 

It  is  this  feature  of  the  Fiacco-McCormick  technique  which  provides 
a  significant  improvement  over  other  convex  programming  techniques. 
Other  approaches,  including  Schmit's  "method  of  alternate  steps"  require 
elaborate  techniques  to  determine  what  action  to  take  when  the  boundary 
is  encountered  during  the  minimization  process. 

Another  consequence  of  the  boundary  repulsion  term  which  should  be 
noted,  however,  is  that  the  search  for  a  feasible  P-function  minimum  must 
be  confined  to  the  region  defined  by  the  inequality  contraints  .     The 
reason  for  this  restriction  is  apparent  from  the  following  example.     Suppose 
that  Gi(X)  =  xx  >  0.    As  xx-»0    from  the  negative  side  [G^X)]"1-  -». 

The  P-function  can  be  minimized  by  direct  analytical  procedures  for 
relatively  trivial  cases  only.    For  other  cases,  some  method  of  successive 
approximations  must  be  used.     One  such  method  is  that  of  gradient 
descent. 
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A.      GRADIENT  DESCENT 

Gradient  methods  are  based  on  a  Taylor's  series  approximation  of 
the  function  to  be  minimized. 

1 .    First-Order  Gradient  Method 

The  first-order  gradient  method  uses  only  the  first-order  terms  of 
the  Taylor's  series  expansion  to  obtain  the  equation, 

X2  =  Xr8  7  P(Xi) 
The  procedure  is  then  to  move  a  distance  proportional  to  9  in  the  direc- 
tion of  the  gradient  vector  evaluated  at  the  point  Xx .    X2    is  the  point  at 
which  the  P-function  is  minimized  along  the  gradient. 
2  .    Second-Order  Gradient  Method 

A  better  approximation  of  the  P-function  is  obtained  by  using  the 
second-order  terms  of  the  Taylor's  series  expansion  as  well  as  the  first 
to  obtain  the  equation 


•  1 


X2  -  Xx  -  Q 


7P(XX) 


dxt  I  Xj 

where  the  matrix  ||  •  ||   is  the  Hessian  matrix  evaluated  at  Xx     and  -1 
denotes  the  matrix  inverse.     The  P-function  is  minimized  along  the  vector 
given  by  the  product  of  the  inverse  of  the  Hessian  and  the  gradient  of  the 
function  evaluated  at  X1# 

If  the  function  to  be  minimized  is  not  strictly  convex,  the  Hessian 
matrix  may  be  positive  definite  only  in  the  vicinity  of  a  local  minimum. 
For  this  reason,  the  procedure  used  began  with  a  search  in  the  direction 
of  positive  6  .    If  the  P-function  did  not  decrease  in  this  direction  before 
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encountering  the  boundary,  the  search  was  continued  in  the  direction  of 
negative  Q  .    This  situation  actually  occurred  in  the  example  when  the 
P-function  was  minimized  at  r=l . 
3  .    Higher-Order  Methods 

It  would  seem  logical  to  consider  the  use  of  higher  order  terms  to 
obtain  an  even  better  approximation  of  the  function  to  be  minimized.    The 
effort  required  to  program  such  a  procedure,  however,  would  be  prohibitive 
for  all  but  trivial  problems.    In  fact,  the  labor  required  to  calculate  the 
second -order  partial  derivatives  for  the  Hessian  matrix  may  prohibit  the 
use  of  the  second-order  gradient  method  for  large  problems  .    Experience 
with  several  small  examples  using  both  techniques  indicates  that  the 
convergence  properties  of  the  second-order  technique  are  far  superior  to 
those  of  the  first-order  technique.     Fiacco  and  McCormick  provide  exam- 
ples of  both  techniques  in  Ref .  3  which  clearly  demonstrate  the  super- 
iority of  the  second-order  method. 

B.      ESTIMATION  OF  THE  FINAL  P-FUNCTION  MINIMUM 

Fiacco  and  McCormick  show  in  Ref.  2  that  estimates  of  the  solution 
of  the  P-function  minimum  for  r=0  can  be  obtained  by  the  use  of  poly- 
nomial approximations  of  x1,x2/  .  .  .  ,x     and  P  as  a  function  of    r£'*     as 
follows  .    Expand  each  component  of  the  X  vector  associated  with  the 
minimum  of  P(X,r)  for  a  given  value  of  r  in  a  power  series  in    r1'2  .    To 
obtain  a  (p-1)  th -order  polynomial  approximation,  drop  the  terms  of  the 
power  series  after  the  pth  term.    If  the  P-function  has  been  minimized  for 
rk ,  k=l ,  .  .  . ,  h  where  h  ^  p ,  form  the  linear  equations : 
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_±±  p-1 

2  2 

x,(r,       )  -  a0  +.  .  .+  a    x(r        )  +  ...+  a        (r        ) 

1    h-p  °  J  1    h-p  p-1    h-p 

•  • 

xt(rh)   =a0  +...+  aj^frh)         +---+ap-i(rh) 

If  it  is  assumed  that  r  is  reduced  by  a  constant  factor  such  that 
rk+1  =  rk  /C  or  rk  =  ra  /Ck_1,    then  the  above  equations  constitute  a  system 
of  p  linear  equations  in  p  unknowns  ,     a0,  .  .  .  , a        .     For  the  polynomial, 


xi  =  a0  +  axr  +  a^r    +  .    .    .  +  a^    ^ 


p-1 

p-1 


xl  at  r=0    is  given  by  a0.     Hence,  if  these  equations  are  solved  for  a0, 
a  (p-1)  th-order  polynomial  approximation  is  obtained  for  the  value  of  xt 
corresponding  to  the  solution  to  the  mathematical  programming  problem. 
Estimates  are  obtained  by  this  procedure  for  x1 ,  .  .  .  ,x     as  well  as  P  at 
r=0.    Because  of  computer  accuracy  limitations,  the  order  of  the  poly- 
nomial estimate  will  be  restricted  to  three  and  will  use  the  last  four 
data  points . 

C.     CRITERIA  FOR  TERMINATION 

As  r  approaches  zero,  the  value  of  the  P-function  minimum  approaches 
the  optimal  solution.     More  rapid  convergence  to  the  optimal  solutiom  may 
be  obtained  by  the  extrapolation  described  above  thereby  permitting  term- 
ination of  the  procedure  with  fewer  iterations.    The  specific  point  at 
which  the  minimization  process  is  to  terminate  depends,  of  course,  upon 
the  accuracy  desired  of  the  final  solution  estimate. 
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Fiacco  and  McCormick  suggest  several  alternate  criteria  in  Ref .  2 
based  on  the  changes  in  the  values  of  the  variables  or  functions  from  one 
P-minimum  to  the  next.    In  the  case  of  the  estimated  value  of  the  P- 
function  minimum,  the  minimization  process  is  terminated  when 

\P*(ri.1)  -  P*  (rt)|<   € 
where  c    is  a  small  positive  number  and  P*(rt)  is  the  polynomial  estimate 
of  the  final  solution  after  the  P-function  is  minimized  for  the  ith  r  value. 
For  the  three-bar  truss  problem,  c  was  chosen  to  be  1x10" 5. 

D.     SUMMARY  OF  SOLUTION  PROCEDURE 

The  following  procedures  are  those  required  to  obtain  an  estimated 
value  for  the  optimal  solution  to  the  mathematical  programming  problem  . 

1.  Form  the  P-function  and  obtain  the  analytical  forms  for  the  gradient 
vector  and  Hessian  matrix. 

2.  Set  p=4  for  a  third-order  polynomial  estimate  of  the  values  of  P  and 
xlt  .  . .  ,x    for  r=0. 

3.  Set  r=l  and  h=l  and  select  a  starting  point,  X0,  such  that  the  in- 
equality constraints  are  satisfied. 

4.  Determine  the  direction  of  steepest  descent  at  the  starting  point. 
Search  along  this  direction  in  the  region  defined  by  the  inequality  con- 
straints until  the  directional  derivative  of  the  P-function  is  less  than  the 
specified  tolerance,  i.e.  near  zero. 

5.  Test  the  gradient  magnitude.    If  less  than  a  specified  tolerance,  go 
to  step  8. 
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6.  Call  the  present  point  the  new  starting  point  and  go  to  step  4. 

7.  If  h    £  p  calculate  the  (p-1)  th-order  polynomial  estimate  of  the  com- 
ponents of  the  solution  at  r=0 . 

8.  If  h    ^  p+1  compare  the  last  two  estimates  of  the  P  function  at  r=0. 
If  within  the  specified  tolerance,  terminate  the  routine.    Otherwise,  set 
r  =  r/C  where  C  is  the  reduction  factor  and  set  h  =  h+1 .    Go  to  step  6. 
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V,     THE  EXAMPLE 
The  example  used  to  illustrate  the  application  of  the  Fiacco- 
McCormick  technique  to  the  design  of  trussed  structures  is  the  problem 
described  previously  with  two  distinct  load  conditions.    This  yields  a 
problem  with  13  variables;  10  equality  constraints,  4  of  which  are  non- 
linear; and  23  inequality  constraints  including  the  three  non-negativity 
constraints  for  the  cross -sectional  areas  of  the  truss  members. 

A.      PROGRAMMING  CONSIDERATIONS 

The  computer  program  was  written  in  Fortran  IV  language  for  the 
IBM  360  computer  at  the  Naval  Postgraduate  School  in  Monterey, 
California.     No  attempt  was  made  to  optimize  the  program  to  reduce 
execution  time.    As  it  is  written,  the  program  required  approximately  2  5 
seconds  for  execution. 

A  simple  bracketing  procedure  was  used  to  locate  the  minimum  of  the 
P-function  along  the  vector  in  the  direction  of  steepest  descent.    The 
directional  derivative  was  evaluated  at  the  current  starting  point.    A    e 
was  then  determined  corresponding  to  a  point  at  which  the  directional 
derivative  was  opposite  in  sign  from  that  at  the  starting  point  and  still 
within  the  region  defined  by  the  inequality  constraints .    The  bracketing 
procedure  was  then  commenced  to  reduce  the  magnitude  of  the  directional 
derivative  below  the  specified  tolerance,  DDTOL. 

Some  experimentation  was  required  to  determine  the  necessary  mag- 
nitude of  DDTOL  to  ensure  convergence  of  the  gradient  magnitude  to  within 


a  tolerance  of  EPSI  =  lxl  0~3.    This  was  determined  to  be  DDTOL  =  lxl  u 
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For  the  three-bar  truss  problem,  a  double-precision  version  of  the 
program  was  required  although  several  two-variable  problems  were 
solved  without  recourse  to  double  precision.    For  small    r    values,  the 
P-function  becomes  progressively  more  peaked  in  the  vicinity  of  the 
minimum  value.    For  the  three-bar  truss  problem  this  peakedness  is  such 
that  the  change  in  ©required  to  obtain  a  bracket  small  enough  to  yield  a 
directional  derivative  within  tolerance  is  below  the  accuracy  of  the 
computer  unless  double  precision  is  used. 

For  programming  purposes,  the  notation  of  the  problem  variables 

W,An  ,A2,A3  ,cj,  ,  ,G^  ,(T31  ,u    ,  ,u    -,  ,cr,  2,000^3  2/u    _,u   n  .  ,  . 

1/    2/    3/    11/    si*    si/    xl/    yl/    12/    22/    32/    x2      y2    was  changed  to 

F,x1,x2,  .  .  .  ,x13    respectively. 

B.      INPUT  DATA 

The  input  data  for  the  problem  was  that  of  Schmit's  Case  1  in  Ref.   1. 

A  cursory  inspection  will  reveal  that  the  values  are  highly  unrealistic 

since  the  actual  range  of  values  of  the  modulus  of  elasticity  is  from  5  to 

30xl06  lbs  ./in.3    for  metals.    They  are,  however,  sufficient  for  purposes 

of  illustrating  the  solution  technique.    The  input  data  is  as  follows. 

Pj.   =  30          Pi    =  1          &  =   135°              E1  =   1 

PS  =  20          pa    =   1          £2  =     90°              E2  =   1 

ax  =  60°         p3    =   1           jS3  =      45°               E3  =   1 

a2  =  180° 

N    =  1 

The  constraints  specified  for  the  problem  are 


25 


A]_  /  A2 ,  A3    ^   0  , 

-15    <  CT1J,a2j,CT3j     <    20,    j  =  l,2, 

and 

-150    <  u    ,,u        <  200,    j=l,2. 
xj      yj 

T 
The  starting  point  for  the  example  was  arbitrarily  chosen  to  be  X    = 

T 
(2,2,2,5,5,5,5,5,5,5,5,5,5)    .    The  reduction  factor  for  r  was  1 0  at 

each  iteration,  i.e.    r1+1=r1/10. 

C.     RESULTS 

Values  obtained  for  minimum  P  as  a  function  of  r  as  well  as  the  third- 
order  estimates  of  the  P-function  minimum  for  r=0  are  shown  in  Table  I. 
The  value  of  the  objective  function,  F,  corresponding  to  the  indicated    r 
value;  the  total  value  of  the  remaining  terms  of  the  P  function,  denoted 
by  P-F;  and  the  values  of  the  equality  constraints,  E.i(X)  are  also  shown. 

The  value  given  under  "Moves  Required  to  Minimize"  is  the  number 
of  times  altew  direction  of  steepest  descent  was  calculated  prior  to 
achieving  a  gradient  magnitude  of  less  than  .001  for  the  current  r  value. 

The  worth  of  the  polynomial  estimate  is  indicated  by  the  fact  that  a 
final  solution  estimate  of  minimum  P  is  obtained  on  the  sixth  iteration 
which  differs  by  less  than  lxl 0"4    from  the  value  obtained  on  the  ninth 
and  final  iteration. 

By  the  ninth  iteration,  the  value  of  the  objective  function,  F,  has 
converged  to  within  7x1 0~!     of  the  final  solution  estimate  while  the  value 
of  P(r9)  differs  by  14xl0~5    from  the  final  estimate. 
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One  consequence  of  not  obtaining  an  absolutely  null  gradient  vector 
is  revealed  by  the  values  of  the  equality  constraints  .    Since  convergence 
is  obtained  only  to  a  gradient  magnitude  of  .001,  the  equality  constraints 
are  not  exactly  satisfied  for  large  r  values.    The  values  of  the  equality 
constraints  are  seen  to  decrease,  however,  at  each  iteration.    In  other 
words,  the  constraints  get  progressively  "tighter"  as  r  -»  0.    The  reason 
for  this  is  the  peakedness  of  the  P-f unction  in  the  vicinity  of  the  minimum 
as  r  gets  smaller  and  smaller.    As  P  becomes  more  peaked,  it  is  necessary 
to  get  nearer  the  minimum  to  satisfy  the  gradient  magnitude  tolerance 
for  convergence.    Thus  the  approximations  become  more  accurate  as 
r  -»   Q.    A    gradient  magnitude  tolerance  of  lxl 0~4  failed  to  yield  a 
change  in  the  value  of  the  constraints  of  sufficient  magnitude  to  change 
the  results  although  it  did  require  a  slightly  higher  execution  time. 

Table  II  shows  the  values  of  the  successive  x^r) ,  i=l  ,  .  .  . ,  13  as 
well  as  the  final  solution  estimates.    The  final  solution  estimates  changed 
by  less  than  2x1 0"4  from  the  sixth  to  the  ninth  iteration. 

The  values  obtained  by  the  method  of  Fiacco  and  McCormick  are 
compared  with  those  obtained  by  Schmit  [1]   using  his  "method  of  alter- 
nate steps"  in  Table  III. 
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Variable 


Ai 


A2 
A3 


all 


°21 


ff31 


U      1 

xl 


u  1 


CT12 


a22 


Q32 


Ux2 


u   0 

y2 


TABLE  III 

Fiacco  & 
McCormick 

1 

.071 

.544 

1 

511 

19 

.877 

20 

.000 

.123 

19 

.755 

20 

.000 

-15 

.000 

5 

.000 

20 

.000 

-35 

.000 

-   5 

.000 

Alternate 
Steps 

1 

.072 

■ 

,544 

■ 

.611 

19 

.863 

19 

.984 

1 

.206 

19 

.743 

-19 

.984 

-14 

.993 

5 

.003 

19 

.996 

-34 

.998 

-   5 

.003 

W  2.922  2.924 

The  results  obtained  by  the  two  methods  are  quite  close  although  the 

minimum  weight  determined  by  the  method  of  Fiacco  and  McCormick  was 

slightly  less.    The  primary  reason  for  the  difference  was  probably  the 

termination  criterion.    Schmit  solved  the  problem  in  1960  on  an  IBM  653 

digital  computer  with  severe  storage  capacity  limitations  so  the  accuracy 

of  the  solution  in  this  study  is  unquestionably  higher.    Advances  in 

hardware  make  any  comparison  of  execution  times  meaningless.     Schmit's 

results  for  o3l   is  erroneous  and  is  probably  a  typographical  error  in  the 

reference. 
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VI.     DISCUSSION 

A.      ALTERNATE  MINIMA 

A  solution  to  the  three-bar  truss  problem  has  been  determined  but 
the  fact  remains  that  it  can  be  verified  to  be  only  a  local  minimum.     Pro- 
cedures exist  by  which  the  "goodness"  of  the  solution  may  be  determined 
but  the  effort  may  be  more  than  that  required  for  solving  the  entire  problem 
by  a  trial  and  error  solution  technique.     One  such  procedure  is  simply  to 
vary  the  starting  point  and  search  for  alternate  minima. 

The  example  problem  was  solved  for  several  other  starting  points 
selected  at  random.    These  included 

X0=   (1. 0,1. 0,1. 0,5.0,5.0,... ,5. 0)T 

T 


.,0.0) 

T 


,0.0) 

T 
.,-5.0) 

..,-io.o)T 


=  (0.5,0.5,0.5,0.0,0.0,. 
=  (0.25,0.25,0.25,0.0,.. 
=  (0.25,0.25,0.25,-5.0,. 
=  (0.25,0.25,0.25,-10.0, 
The  only  change  in  value  noted  was  a  difference  in  the  number  of 
moves  required  to  minimize  the  P-function  for  r=l . 

B.      TWO-BAR  TRUSS 

If  the  center  bar  of  the  truss  is  removed,  the  structure  becomes  a 
statically  determinant  two-bar  truss  in  which  the  internal  forces  are 
uniquely  determined.    The  minimum  weight  of  such  a  structure  subject  to 
the  constraints  of  the  example  is  3.04905.    Thus,  by  adding  a  third 
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member,  the  total  weight  of  the  structure  is  reduced.    It  should  be  noted, 
however,  that  the  behavior  of  the  trusses    will  not  be  the  same. 

Switsky  [4]  shows  that  the  deflection  of  the  two-bar  truss  will  be 
less  in  this  case  than  that  of  the  three-bar  truss.    If  the  displacements 
of  both  trusses  are  required  to  be  equal,      given  the  same  stress  con- 
straints, the  statically  determinant  two-bar  truss  will  be  the  lighter  of 
the  two. 

C.     NON-ANALYTICAL  CONSTRAINTS 

Structural  problems  often  do  not  allow  the  nice  problem  statement 
form  of  the  three-bar  truss  problem.    The  major  difficulty  arises  from  not 
being  able  to  express  some  of  the  constraints  in  the  completely  analytical 
form  required  for  the  Fiacco-McCormick  technique.    Many  problems 
require  the  use  of  empirically  determined  alignment  charts  or  graphs  to 
calculate  changes  in  structural  behavior  as  a  result  of  the  variation  of 
design  variables.    For  such  problems,  some  variation  of  Schmit's 
"method  of  alternate  steps"  may  be  well  suited.    An  example  of  such  a 
problem  is  the  determination  of  the  minimum  weight  design  of  an  integrally 
stiffened  panel  loaded  axially  under  compression.    Rosenbaum  [5]     des- 
cribes the  problem  and  outlines  a  solution  procedure  based  on  the  method 
of  alternate  steps. 
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VII.     RECOMMENDATIONS  FOR  FURTHER  STUDY 

Since  only  three  design  variables  of  the  three-bar  truss  are  allowed 
to  vary,  the  example  is  actually  a  sub-optimization  problem.    This  fact 
suggests  that  several  possible  extensions  of  the  problem  merit  further 
consideration. 

Allowing  more  design  variables  such  as  the  length  of  truss  members 
to  vary  makes  the  problem  even  more  nonconvex  but  the  Fiacco-McCormick 
procedure  may  still  be  applied.    The  assumption  of  the  availability  of  a 
continuous  spectrum  of  materials  would  permit  the  determination  of 
optimal  material  properties  of  the  structure. 

The  added  restriction  of  a  discrete  choice  of  materials  or  material 
sizes  available  is  a  further  realistic  extension  of  the  problem  although  a 
different  solution  procedure  is  required. 

Finally,  the  development  of  a  systematic  procedure  for  finding  alter- 
nate minima  is  essential  if  one  is  to  accept  any  result  with  confidence. 
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COMPUTER  PROGRAM 

C    MAIN    PROGRAM 
r 

IMPLICIT    FEAL*B     (A-H,0-Z) 

REAL*B  KllfK12fKl3tK21 t K22  ,K2  3 , K1A ,K1 5 ,K2  5 ,K31 ,K42 ,K5  3 
$  ,  K  34  ,  K  35  t  K  W  ,  K  AS  t  K  54 ,  K  55 .  K 24 

DIMENSION  X(13),XX(13,12),XEST(13,12),X2U3),SPR(13,13 
$),COL(  lfc<5)  .OELPt  13)fXDLO(  13)  ,LWK(  13)  ,MWK(  1  3)  ,C  ( 1 3  )  ,H  ( 1 
$C  ,12) 

DIMENSION    RR(12),P<12),F<12) , SIGMA (12 ),PFST(12),IA(12) 

DIMENSION     XY(4),RY(4) ,PY(4) 

F QUI  VALENCE ( SPR( 1  ) ,COL(l ) ) 

COMMON    Cl,C2»C3,Kll,K12,K13,K21,K2  2,K2  3fK14,K24,K15,K2 
$5,K31 ,  K42, K53,K34tK3  5, K44, K45,K54,K55 
C    QDTOL    -    TOLERANCE    FOR    DIRECTIONAL    DERIVATIVE. 
C    EPSI    -    TOLEFANCE    FOR    GRADIENT    MAGNITUDE 

DATA    DDTOL/.C?r"C01/,EPSI/.Cm/tDELTA/l./ 
C     STARTING    X    VALUES    FOR    SEARCH. 

DATA    X/3*-2.:,ir*5.G7 
C    NV    -    NUMBER    OF    VARIABLES. 

NV=13 

RF  (1  )=1. 
C    IR    -    NUMBER    OF    !*    VALUES    FOR    WHICH    P    FUNCTION     IS    TO    BE 
C    MINIMIZED. 

IP=9 
C    J    -    NUMBER    OF    P    VALUES    FOR    WHICH    P    HAS    BEEN    MINIMIZED. 

J  =  ^ 
5    J=J+1 
C    RR (J)     -    JTH    P    VALUE. 

R=RR( J) 
C 

C  BEGIN  SECOND  ORDER  GRADIENT  SEARCH  FOR  MINIMUM  OF  P  FOR 
C  CURRENT  R  VALUE. 
C 
C  DDOLD  -  PREVIOUS  DIRECTIONAL  DERIVATIVE. 

DDOLD=O.C 
C  GOLD  -  PREVIOUS  GRADIENT  MAGNITUDE. 

GCLD=1.C 
C  XOLD  -  X  VALUFS  CORRESPONDING  TO  PREVIOUS  GRADIENT 
C  MAGNITUDE. 

DO  11  1=1, NV 

XOLD( I )=f.C 
11  CONTINUE 
C  II  -  COUNT  OF  MOVES  IN  SEARCH  FOR  MINIMUM  OF  P  FOR  CURRFNT 
C  R  VALUE. 

11=0 
C  EVALUATE  HESSIAN  MATRIX  AT  CURRFNT  STARTING  POINT. 
C  SPR  -  HESSIAN  MATRIX. 

13  CALL  MSOPAPU,  R,SPR,NV) 
C  EVALUATE  GRADIENT  VECTOR  AT  CURRENT  STARTING  POINT. 
C  DELP  -  GRADIENT  OF  P. 

CALL  GRAD( X,R,DELP) 
C  INVERT  HFSSIAN  MATRIX. 

CALL  MINV(COL,NV,D,LWK,MWK) 
C  STOP  IF  HESSIAN  MATRIX  IS  SINGULAR.  PICK  NEW  STARTING 
C  «»OINT. 

IFID.EQ.r ,C )  WRITE(6t 17) 
17  FtRMAT(«*DET  OF  MATRIX  OF  2ND  ORDER  PARTIALS  =  ?*•) 

IF(D.EQ.r.C)  STOP 
C  PREMULTIPLY  INVERSE  OF  HESSIAN  BY  GRADIENT  VECTOR. 

CALL  MTXMUL(SPR,DELP,NV»C) 
C  GMAG  -  MAGNITUDE  OF  GRADIENT. 

GMAG=C.C 

DO  2  3  I=1,NV 

GMAG=GMAG+(DELP(  I )*OELP( I )  ) 
23    CONTINUE 

GMAG=DSQPT(GMAG) 
C    TEST    MAGNITUDE    OF    GRADIENT.     IF    BELOW    TOLERANCE,    P     IS 
C    MINIMIZED    FOR    CURRENT    R. 

IF(GMAG.GT.EPSI )    GO    TO    28 
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Y=DD/DAB 
(TRY+TES 

OPPOSIT 
M  IN  DIR 
(PV2.LT. 
EP=-2. 

TO  52 


IONAL  DERIVATIVE.  KEEP  SIGN.  EVALUATE  P 

ART  POINT. 

BS(DD) 

(X,R) 

VALUE  IS  POSITIVE. 


ION  OF  STEEPEST  DESCENT  TO  GET  X2. 

E( X2,X,NV,C,T) 

ASIBILITY    WITH    RESPECT    TO    INEQUALITY 

(X2,NV, RESULT) 

NE,  X2  IS  FEASIBLE.  IF  X2  IS  NOT  FEASIBLE, 

ALUE. 

EQ.1.0)  GO  TO  58 

2. 


ER  STEP. 

(  X  2  R  ) 

TIONAL  DERIVATIVE  AT  X2. 

X2,R,DELP) 

NV,DELP.C) 

UM  AT  X2. 

).LE. DDTOL)  GO  TO  82 

IRECTIONAL  DERIVATIVE  AT  X2.  IF 

TARTING  POINT,  BRACKET  HAS  BEEN 

S(DD> 

T)  65,63,65 

E  IN  SIGN  BUT  P  INCREASES,  SEARCH  FOP 

ECTION  OF  NEGATIVE  THETA  FROM  START  POINT. 

PV1)  GO  TO  52 


OPPOSITE 
ESTABLISHED. 
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C    INITIAL    BRACKET    IS    ESTABLISHED.    REDUCE    STEP    SIZE    BY    1/2. 
C    REVERSE    DIRECTION    OF    SEARCH    WITH    REDUCFD    STEP    SIZE. 
68    TEST=TRY 
STEP=T 

70  STEP=-STFP/2. 

71  T=T+STEP 

CALL    XVALUE(X2,X,NV,C,T) 

CALL    GRADI X2,R,DELP) 

DD=DIROER(NV,DELP,C) 
C    TEST    FOP    MINIMUM    AFTER    EACH    STEP. 

IF(DABS(PD).LF.DDTOL)     GO    TO    82 
C    TERMINATE    SEARCH    IF    DIRECTIONAL    DERIVATIVE    REPEATS. 

IF(DD.EQ.DDOLD)    GO    TO    82 
C    SAVE    CURRENT    VALUE    OF    DD. 

DDCLD=DD 
C    TEST    FOR    ESTABLISHMENT    OF    REDUCED    BRACKET. 

TPY=DD/DABS(DD) 

IFITRY+TEST)     71,80,71 
80    TEST=TRY 

GO    TO    7C 
C    CALCULATE    VALUE    OF    NEW    STARTING    POINT.    BEGIN    SEARCH    IN 
C    NEW    DIRFCTION    OF    STEEPEST    DESCENT. 
82    CALL    XVAt UE( X , X ,NV ,C , T ) 

11=11+1 

GO    TO    13 
85    I A  ( J  )=  1 1 

P( J)=PSTAR 
C    EVALUATE    ORIGINAL    OBJECTIVE    FUNCTION. 

F< J)  =  X  (1 )*C1+X<2)*C2+X(3)*C3 

SIGMA( J)=P<J)-F{ J) 
C    EVALUATE    EQUALITY    CONSTRAINTS. 

H( 1, J)=X(1 )*X(4)*KL1+X(2)*X(5)*K12+X(3)*X< 6)*K13-K14 

H<  2,  J)  =X(  1  )*X<4)*K21  +  X(2)*X(5)*K22+xm*X<6)*K23-K24 

H(?,J)=K31*X(4)+K34*X17)+K35*X<8) 

H(4,  J)=K42*X(5)+K44*X(7)-»-K45*X(8) 

H(  5,  J)=K53*X(  6)+K54*X(7)+K55*X(  8) 

H<6»J)=X(1 )*X(9)*K1H-X(21*X(10)*K12+X( 3) *X < 1 1 ) *K 13-K15 

H(7, J) =X(1 )*X(9)*K2H-X(2)*X(1C)*K22  +  X(3)*X (ll)*K23-K2  5 

H(  8,  J)=K31*X(9  )+K34*X<  12)  +K35*X  (  1  3  ) 

H<9, J)=K42*X(10)+KVf*X(12)+K45*X(  13) 

H< 1C,J)=K5  3*X( 11)+K54*X(12)*K55*X( 13) 

DO    9  2    K=l, NV 

XX(K, J)=X( K) 
92    CONTINUE 

IFU.LT.4)     GO    TO    108 
C 

C    BEGIN    THIRD    ORDER    POLYNOMIAL    CURVE    FIT    TO    ESTIMATE    P    AND    X 
C    VALUES    FOP    R=C  .    PY,RY,AND    XX    CONTAIN    LAST    FOUR    VALUES    OF 
C    P,R,AND    X     IN    FCKM    REQUIRED    FOR    SUBROUTINE    POLY. 


C 


11=0 

J3=J-3 

DO    100    K=J3,J 

11=11*1 

PY(I1)=P(K) 

RY(I1)=RR( K) 
130    CONTINUE 

CALL    PCLY(PY,RY,4, AO) 
PEST    -    ESTIMATE    OF    P    AT    R=0. 

PEST(J)=AO 

DO    105    1=1, NV 

11=0 

J3=J-3 

DO    102    K=J3,J 

11  =  1  1+1 

XY(I1)=XX( I,K) 
102    CONTINUE 

CALL    POLY( XY,RY,4,A0) 
XEST    -    ESTIMATE    OF    X    AT    R=0. 

XESTd  ,J)  =  AO 
105   CONTINUE 

IF(J.EQ.IR)    GO    TO    122 
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108    J1=J*1 
C    REDUCE    P    BY    ARBITRARY    CONSTANT. 
RR(J1)=RR(  JI/10. 
GO    TO    5 
122    WRITE(6,123) 

12  3    FORMAT ( • 1 • 44  X • MOVES • / 4  3X ' REQU I  RED ■ /46X «TO • /23X«L I NE • 8X 
$»R»7X»MIM  MIZE • 4X» P« 8X • F • 7X» P-F' 8X«H1 « 8X • H2» 8X »H3» ) 
DO    133    J=l tIP 

WPITE(6,128)J,RR(J),IA(J),P(J),F(J ), SIGMA (J), (H( I, Jit  I 
$=1,3) 
128    FORMAT ( • C •  2?X,  13 , 3 X , F 1 2.9 , 3X , I  3 t 3X ,F8. 5 , F9 .5 , F 10 . 5 , 
$F1C.5,F1C. 5,F10.5) 
IF(J.LT.A)    GO    TO    133 
WPITE(6,132>     PEST(J) 

132  FORMATS     »29X'3RD    ORDER    ESTIMATE       »F8.5) 

133  CONTINUE 
WRITE(6,135) 

135  FORMAT (• 1 • 23 X • L I NE • 7X» R' 1AX* H4« 8X • H5 • 8X» H6 »3X • H7 » 8X »H8 
$8X,H9«8X»H10»  ) 

WRITE(6,138)(J,RR(J),(H(I,J),l=4,K)tJ  =  lfIR) 
138  FORMAT  PC  23X,I3,F14.9,7F10.6) 

WRITE(6,14C) 
140  FORMAT  (•  l»23X«LINE,7X»R,llX,Xl»8X,X;>,aX,X3,3XfX4«8X,X5 
$8X'Xb9 8X»X7' ) 

DO  148  J  =  l  ,IR 

WRITE (6, 144) J,RR( J) , (XX( I, J) , 1=1,7) 
144  FORMAT <«C • 23X , I  3, F 14. 9, 7F 10. 5 ) 

IFIJ.LT.4)  GO  TO  148 

WRITEI6.147)     ( XESTI I, J) ,1=1,7) 

147  FORMAT!*     »,28X»3RD    ORDER    ES t • F9. 5 , 7F 10 . 5 ) 

148  CONTINUE 
WRITE(6,15C) 

150  FORMAT (• 1«  28X» LINE* 7X«R' 1 IX' X8»  3X,X9«8X,X10»7X,X11«7X 
$,X12,7X»X13' ) 

DO  158  J=1,IR 

WPITE(6,154)J,RR(J),<XX(I,J),I=8,13) 
154  FCiRM  AT  (•("■•  23X,I3,F14.9,6F10.5) 

IFU.LT.4)  GO  TO  158 

WRITE (6, 15 7)  (XEST(  I, J), 1  =  8, 13) 

157  FORMAT!'  ',33X'3RD  ORDER  EST • F9. 5 , 6F 10 .5  I 

158  CONTINUE 
STOP 
END 

C 
C 
C  SUBPROGRAMS 

SUBROUTINE  POLY ( X , RR, K , AO ) 
C  THIS  SUBROUTINE  SOLVES  THF  LINEAR  EQUATIONS  FOR  THE 
C  POLYNOMIAL  FSTIMATE  OF  THE  VARIABLE  CONTAINED  IN 
C  ARGUMENT  X.   AO  IS  THE  ESTIMATED  VALUE  OF  THE  VARIABLE 
C  AT  R-r . 

IMPLICIT  R EAL*8(A-H,0-Z) 

DIMENSION  X(4),RR(4),A(4,4),L(16),M(16),C(16) 

EQUIVALENCE!  A(  1)  ,C(1) ) 

DO  1  J=1,K 

DO  1  1=1, K 

A( I, J) =<DSQ*T(RR( I ) ) )**( J-l 5 

1  CONTINUE 

CALL  MINV(C,K,D,L,M) 

AO=C  .0 

DC  2  J1=1,K 

A0  =  A0-MA<1  ,J1)  )*X(J1) 

2  CONTINUE 
RFTURN 
END 

C 
C 

DOUBLE  PRECISION  FUNCTION  DIRDER ( NV, DELP ,C  ) 
C  THIS  SUBPROGRAM  CALCULATES  THE  VAL'»E  OF  THE  DIRECTIONAL 
C  DERIVATIVE. 
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IMPLICIT    REAL*8(A-H,0-Z) 
RFAL*8    DFLP(NV) ,C(NV) ,MAGSQ 
MAGSQ=O.C 
DO    1     1=1, NV 
MAGSQ=MAGSQ-HC(I  )*C(  I)  ) 

1  CONTINUE 

IF     (MAGSQ.EQ.O.n)    GO    TO    3 

UV=1/DSQRT<MAGSQ) 

DIRDER=C.O 

DO    2    1=1, NV 

DIRDER=DIRDER+(DELP(I )*C(I )*UV) 

2  CONTINUE 
RETURN 

3  DIRDER  =C.C 
RETURN 

END 

SUBROUTINE    MTXMUL ( SPR ,DELP, NV, C ) 
C    THIS    SUBROUTINE    PREMULTIPL I ES    THE    VFCTOR    DELP    BY    THE 
C    MATRIX    SPR. 

IMPLICIT    PEAL*8(A-H,0-Z) 

DIMENSION    SPR(NV,NV),DELP(NV),C(NV) 

DO    1     1=1, NV 

C(  I)=C.C 

DO    1    J=1?NV 

C(  I)=C(  I)-MSPR(I,J)*DELP(  J)  ) 
1    CONTINUE 

RETURN 

END 

SUBROUTINE    XVALUE     ( X2 ,X1 , J,C, T  ) 
C    THIS    SUBROUTINE    CALCULATES    THE    POINT    X2    A    DISTANCE 
C    PROPORTIONAL    TO    T    FROM    XI     IN    THE    DIRECTION    OF 
C    STEEPEST    DESCENT. 

IMPLICIT    REAL*8(A-H,0-Z) 

DIMENSION    X2( J),Xl(J) ,C(J,  1) 

DO    1    I     =    1,J 

X2(I )    =    XI (I )-(T*C<  1,1)) 
1    CONTINUE 

RETURN 

END 

SUBROUTINE    CHECK ( X ,NOVAR, RESULT  ) 
C    THIS    SUBROUTINE    CONTAINS    THE     INEQUALITY    CONSTRAINTS. 
C    IF    THE    POINT    X    SATISFIES    THESE    CONSTRAINTS,    RESULT=+1. 

IMPLICIT    REAL*8(A-H,0-Z) 

REAL*8  X(13),LA,L5,L6,L7,L8,L9,L1C,L11,L12,L13,K11,K12 
$, K 1 3, K 2 1, K 22, K2 3, K 14, K 1 5, K25,K 3 1,K 42, K 53, K 34, K35,K4A, 
$K45,K54,K55,K24 

COMMON  Cl,C2,C3,Kil,K12,K13,K21,K22,K23,K14,K24,K15, 
$K2  5,K31,K4  2,K5  3,K34,K3  5,K44,K45,K54,K5  5,U4,U5,U6,U7,U8 
$,L4,L5,L6,L7,L8 

EQUIVALENCE  (L4,L9),(L5,L1C),(L6,L11),(L7,L12),(L8,L13 
$) ,(U4,U9), (U5,U10) ,(U6,U11) ,(U7,U12), (Ufl,U13) 

IF(X( 1).LE .3.0)  GO  TO  1 

IF(X(2).LE.Q.O)  GO  TO  1 

IFIxm.LE.O.O)  GO  TO  1 

IF(X(4).LE.L4.0R.X<4) .GE.U4)  GO  TO  1 

IF(X(5).LE.L5.0R.X<5) .GE.U5)  GO  TO  1 

IF(X(6).LE.L6.0R.X(6) .GE.U6)  GO  TO  1 

IF(X(7).LE.L7.0R.X(7).GE.U7)  GO  TO  1 

IF<X(8).LE.L8.0R.X(8) .GE.U8)  GO  TO  1 

IF (X(<M.LE.L9.0R.X<9) .GE.U9)  GO  TO  1 

IF(X<1C)  .LE.L10.OR.X(  lO.GE.UlD)  GO  TO  1 

IF(X(11) .LE.Lll.OR.X(ll) .GE.U11)  GOTO  1 

IF(X(12) .LE.L12.0R.X1 12).GE.U12)  GO  TO  1 

IF(X(13).LE.L13.0R.X( 13J.GE.U13)  GO  TO  1 

RESULT=1.0 

RETURN 
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1    RESULT=-1.0 

RETURN 

END 

DOUBLE    PRECISION    FUNCTION    PVALUE(X.R) 
C    THIS    SUBPROGRAM    CONTAINS    THE    ANALYTICAL    FORM    OF    THE 
C    P-FUNCTION.       THE    VALUE    OF    P    AT    X    I S    CALCULATED. 

IMPLICIT    REAL*8(A-H»0-Z) 

REAL*?    X(13),LA,L5,L6tL7,L8,L9tL10tLlltL12,Ll3tKll,K12 
$,K13,K21,K2?iK23,K14,K15,K25,K31,K42,K53,K34,K35,K44, 
$K^5fK54,K55fK2A 

COMMON    C1,C2,C3,K11,K12,K13,K21,K22,K23,K14,K24,K15, 
$K2  5,K31,K4  2,K5  3,K34,K35,K44,K45,K54,K55tU4,U5,U6tU7,U8 
$tLA,L5,L6, L7.L8 

EQUIVALENCE    (L4,L9)f<L5,L10),(L6,Lll),(L7,Ll2),(L8,L13 
$) .  (U4,U9), (U5tU10l , (U6fUll) t (U7,U12), (U8,U13) 

Hl=X(I)*X(4)*Kll+X(2)*X(5)*K12+X(3 )*X (6 ) *K 13-K14 

H2  =  X(  1  )*X<  4)*K21*X(2)*X(5)*K22*X(3)*X(6)*K23-K24 

H3  =  K31*X(4)+K34*X(7H-K3  5*X(8) 

H4  =  K42*X(5)  +  K44*X(7)-»-K45*X(  8) 

H5=K53*X<6  )+K5  4*X(  7)  ♦  «  55*X(  8  ) 

H6  =  X(1  )*X(9»*K11«-X(2I*X(10>*K12+X(3)*X(11)*K13-K15 

H7=X( 1)*X(9)*K21+X(2)*X( 10)*K2  2+X( 3 ) *X ( 1 1 ) *K23-K2  5 

H8=K31*X(9 )+K34*X( 12)+K35*X( 13) 

H9=K42*X( 1C)+K44*X(12)+K45*X(13) 

H1C=K5  3*X(  11)«-K54*X( 1 2 ) ♦K55*X< 1 3) 

Gl=l./X(l ) 

G2=l./X(2) 

G3=l./X(3) 

G4=1./(X(4)-L4) 

G5=1./(X(4)-U4> 

G6=1./(X(5)-L5) 

G7=1./(X(5)-U5) 

G8=1./(X(6)-L6) 

G9=l./ (X(6 )-U6 ) 

G10=1./(X( 7J-L7) 

G11=1./(X( 7)-U7) 

G12=1./(X(8)-L8) 

G13=1./<X< 8J-U8) 

G14=1./(X(9)-L9) 

G15=1./(X(9)-U9) 

G16=1./(X( 1 JI-L10) 

G17=1./<X( 10I-U10) 

G18=l.  /(X(  1D-L11I 

G19=1./(X(  1D-U11) 

G20=l./(X< 12)-L12) 

G21=1./(X( 12I-U12) 

G22=1./<X( 13J-L13) 

G23=1./(X( 13)-U13) 

Y1=1./DSQPT(R) 

EQ=Hl**2-«-H2**2-»-H3**2  +  H4**2  +  H5**2  +  H6**2+H7**2-»-H8**2-«- 
$H9**2+H10**2 

UN  =  GH-G2+G3+G4-G5*-G6-G7+G8-G9+Gin-Gll  +  G12-G13  +  G14-G15  + 
$G16-G17+G18-G19+G20-G21+G2  2-G2  3 

PVALUE=C1*X( 1 UC2*X(2)+C3*X(3)+Y1*EQ+R*UN 

RETURN 

END 

SUBROUTINE    GRAD(X,R,D) 
C    THIS    SUBPROGRAM    CONTAINS    THE    ANALYTICAL    FORM    OF    THE 
C    GRADIENT    VECTOR.       THE    VECTOR    IS    FVALUATED    AT    X. 

IMPLICIT    REAL*8(A-H,0-Z» 

REAL*8    DI13) 

PEAL*8    X(13),L4,L5,L6,L7,L8,L9,L1C,L11,L12,L13,K11,K12 
$f K13,K21tK22tK23tK14,Kl5,K25fK31,K42fK53,K34,K35,K44, 
$K45,K54,K55,K24 

COMMON    CltC2tC3,Kll,K12,K13,K21,K22,K23,K14,K24,K15, 
$K2  5,K31,K4  2tK5  3tK3AfK3  5tK44,K4  5,K54,K55fU4,U5tU6,U7fU8 
$,L4,L5,L6,L7,L8 

EQUIVALENCE    (L4,L9) . ( L5 ,L1C) . ( L6 ,L11 » , ( L7, L12 » , ( L8 ,L13 
i)i  (U4,U9),  (U5,U10),(U6tUll),  (U7.U12)  ,(U8,UU) 

DIMENSION    RO(3) , E ( 3 ) , BA( 3 ) , AA( 2 ) , P (2 ) 
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H1=X(1)*X(4)*K1H-X(2)*X(5)*K12*X(3)*X(6)*K13-K14 

H2=X(1  )*X14)*K2H-X(2)*X<5)*K22*X<  3  )*X(  6)  *K  23-K24 

H3=K31*X(4)*K34*X(7)+K35*X(8) 

H^=KA2*X(  5  )+K4  4*X(  7)  +K45*X(  3) 

H5=K53*X(6)+K54*X(7)*K55*X(8) 

H6=X(1)*X< 9)*K11+X(2)*X(10)*K12+X(3)*X(11)*K13-K15 

H7=X(1)*X(9)*K2H-X<2)*X< 1C)*K22+X( 3 )*X< 1 1 ) *K23-K2  5 

H8=K31*X(9)*K34*X(12)-«-K35*X<13) 

H9=K42*X( lC)+K44*X(12)-»-K45*X( 13) 

H10=K5  3*X(  ll)-»-K54*X<  12H-K55*X<  13) 

Gl=l./X(l) 

G2=l./X(2) 

G3=l./X(3) 

G4=1./(X(4)-L4) 

G5=1./(X(4)-U4) 

G6=1./(X(5 )-L5) 

G7=1./(X<5)-U5) 

G8=1./(X(6)-L6) 

G9=1./(X(6)-U6) 

G10=1./(X(  7)-L7) 

G11=1./(X(7)-U7) 

G12=1./(X(8)-L8) 

G13=1./(X(  8)-U8) 

G14=1./(X(9)-L9) 

G15=1./(X( 9)-U9) 

G16=1./<X<  10J-L10) 

G17=1./(X(  1O-U10) 

G18=l. /<X<  11 )-Lll ) 

G19=1./(X<  1D-U1D 

G2C=1./(X( 12)-L12) 

G21=l.  /(X(  12)-U12) 

G22=1./(X( 13)-L13) 

G23=l. /<X< 13J-U13) 

Y2=2./DSQPT(R) 

D(  l)=Cl+Y2*(X(4)*(Kil*HH-K21*H2)  +  X(9)*(Kll*H6-»-K2l*H7) ) 
4— R*(  Gl  ♦♦?  1 

D(2)=C2+Y2*( X( 5)*(Kl2*HH-K22*H2)*X(10J*(Kl2*H6+K22*H7i 
$)-R*(G2**2 ) 

D<  3)=C3-»-Y2*(X(6)*(K13*HH-K23*H2H-X<li  )  *<  Kl  3*H6*K2  3*H7  ) 
*)— P*(G3**? ) 

D(4}=Y2*(X<l)*(Kli*Hl+K21*H2)+K3l*H3)-R*(G4**2-G5**2) 

0(  5)=Y2*(X(2)*(K12*H1+K22*H2)+K42*H4)-R*(G6**2-G7**2) 

0(6)=Y2*(X (3)*<K13*H1+K23*H2)+K53*H5)-R*(G8**2-G9**2) 

D( 7)=Y2*(K34*H3+K44*H4+K54*H5)-R*(G10**2-G11**2) 

D<8)=Y2*(K35*H3+K45*H4+K55*H5)-R*<G12**2-G13**2) 

0<9)=Y2*(X  <1)*(K11*H6*K21*H7)+K31*H8)-R*(G14**2-G15**2 
$) 

D(  1C  )  =  Y2*(  X<  2)*(K12*H6+K22*H7)+K42*H9)-R*(G16**2- 
$G17**2) 

0( 11)=Y2*< X(3)*(K13*H6*-K23*H7)«-K53*H10)-R*(G18**2- 
$G 1 9**2  ) 

D(12)=Y2*(K34*H8+K44*H9+K54*H10)-R*(G20**2-G21**2) 

D( 13)=Y2*< K3  5*H8+K45*H9+K55*H10)-R*<G2  2**2-G23**2) 

RETURN 

END 

SUBROUTINE    MSOPAR ( X.R , S.NV I 
C    THIS    SUBROUTINE    CONTAINS    THE    ANALYTICAL    FORM    OF    THE 
C    HESSIAN    MATRIX.       THE    MATRIX    IS    EVALUATED    AT    X. 
C    PROBLEM    CONSTANTS    ARE    READ    INTO    COMMON    IN    THIS    ROUTINE. 

IMPLICIT    REAL*8(A-H,0-Z) 

REAL*8    S(13, 13) 

REAL*8    N 

REAL*8    X(13) ,L4,L5.L6,L7,L8,L91L10,L111L12,L13,K11,K12 
$tK13,K21,K22,K23,Kl4,K15,K25,K3l,K42,K53,K34,K35,K44, 
$K45tK54,K55.K24 

COMMON    Cl,C2tC3,KlltK12tK13,K21»K22tK23tK14,K24,K15, 
$K25tK31,K42tK5  3tK34,K35tK44,K45,K54fK55tU4,U5,U6,U7,U8 
$,L4fL5,L6,L7,L8 

EQUIVALENCE    (L4,L9)t(L5,L10)t<L6,Lll).<L7fL12)i<L8,L13 
S) ,<U4,U9), (U5.U10) , (U6,U11), (U7,U12), (U8,U13) 
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01  HENS  ION 

DATA    RO(l) 

DATA    L/C/ 

L  =  L+1 

IF(L.GT.l) 

U4=20.0 

05=20.0 

06=20.0 

07=200. 

U8=200. 

L4=-15. 

L5=-15. 

U=-15. 

L7=-2CC. 

L8=-200. 

PI=3. 14159 

BA(1)=.75* 

BA(2)=.5C* 

BA<3)=.25* 

AA(l)=PI/3 

AA(2)=PI 

P<1)=3C. 

P(2)=2C. 

C1=N*R0< 1) 

C2=N*P0<2) 

C3=N*R0<3) 

K11=DC0S(B 

K12=DC0S(B 

K13=DC0S(B 

K21=DS  IN(B 

K22=DSIN(B 

K23=DSIN(B 

K14=-P(l)* 

K24=P( 1)*D 

K15=-P(2)* 

K25=P(2)*D 

K31=N/(E(1 

K42=N/(E(2 

K53=N/tE(3 

K34=K11 

K35=K21 

K44=K12 

K45=K22 

K54=K13 

K55=K23 

Y2=2./DS0R 

H1=X(1)*X( 

H2  =  X(1  )*X( 

H3=K31*X(4 

H4=K42*X( 5 

H5=K53*X<6 

H6=X(1 )*X< 

H7  =  X(1  )*X( 

H8=K31*X(9 

H9=K42*X( 1 

H10=K53*X( 

Gl=l./X(l) 

G2=l./X(2) 

G3=l./X<3) 

G4=1./(X<4 

G5=1./(X(4 

G6=1./(X(5 

G7=1./(X(5 

G8=1./<X<6 

G9=1./(X(6 

G10=1./(X( 

G11=1./(X( 


R0(3) ,E(3) ,BA(3) .AA(2),P<2) 

iRO(2)tRO(3),E(l)iF(2litnitN/7*l.V 


GO    TO    2 


26536 
PI 
PI 
PI 


/0SIN(BA(1 ) ) 
/0SIN(BA(2) ) 
/0SIN(BA(3) ) 
A(  1)  ) 
A(2)  ) 
A<3)  ) 
A(l)  ) 
A(2)  ) 
A(3)  ) 

DCOS(AA(  1)  ) 
SIN<  AAU  )) 
DC0S(AA(2) ) 
SIN( AA(2)) 
)*DSIN(BA<1) ) ) 
)*DSIN(BA(2)  )  ) 
)*DSIN(BA<3) ) ) 


G12=l 
G13=l 
G14=l 
G15=l 


/(X( 
/(X( 
/(X( 
/(X( 


G16=1./(X( 


T(R) 

4)*K11+X(2)*X(5)*K12*X(3)*X<6)*K13-K14 

4)*K2H-X(2)*X(5)*K22  +  X<3)*X(6)*K2  3-K24 

)+K34*X(7)+K35*X(3) 

)+K44*X(  7)+K45*X(  3) 

)+K54*X(7)+K55*X(8) 

9>*Kll+X(2)*X(10)*Ki2+X(3)*X(H)*K13-K15 

9)*K21+X(2)*X(  1C)*K2  2*X(3>*X(11)*K23-K25 

)+K34*X< 12)+K35*X< 13) 

C)+K44*X<12)*K45*X< 13) 

11)*K54*X< 12)+K55*X(13) 


)-L4) 

)-04) 

)-L5) 

)-U5) 

)-L6) 

)-U6) 

7)-L7) 

7)-U7) 

8)-L8) 

8)-U8) 

9)-L9) 

9)-U9) 

10J-L10) 
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G17=l. 
G18=l. 
G19=l. 
G20=l. 
G21=l. 
G22=l. 
G23=l. 
S(l,  1) 
S< 1,2) 
S(l,3) 
S< 1,4) 
S(l,5) 
S( 1,6) 
S(  1,7) 
S(l,8) 
S( 1,9) 
S(  1,1C 
S  (  1 , 1 1 
S(  1,12 
S( 1,13 
S(2,2) 

$G2**3 
S(2,3) 
S(2,4> 
S(2,5) 
S(2,6) 
S(2,7) 
S(2,8) 
S<2,9) 
S( 2,10 
S(2,ll 
S(2,12 
S( 2,13 
S(3,3) 

$G3**3 
S<3,4) 
S(3,5) 
S( 3,6) 
S(3,7) 
S(3,8) 
S<  3,9) 
S<3,1C 
S(3,ll 
S(3,12 
S(3,13 
S(4,4) 

$(G4**3 
S<4,5) 
S(4,6) 
S<4,7) 
S(4,8) 
S(4,<U 

S(4,10 
S<4,11 
S<4,12 
S(4,13 
S( 5,5) 

$(G6**3 
S(5,6) 
S(5,7) 
S(5,8) 
S(5,9) 
S(5,10 
S(5,ll 
S(5,12 
S(5,13 
S(6,6) 

$(G8**3 
S(6,7) 
S(6,8) 
S(6,9) 
S(6,1C 


12 
12 
13 
13 


/<X<  1) 

/(X(ll 

/(X(  11 

/<X< 

/(X( 

/(X( 

/(X( 

=Y2*(X 

=Y2*(X 

=Y2*  (X 

=Y2*<K 

=Y2*X( 

=Y2*X< 

=0.0 

=0.0 

=Y2*(K 

)  =  Y2*( 

)=Y2*( 

)=C.O 

)  =  C.C 

=Y2*  (X 


)-U10) 

)-Lll) 

)-Ull) 

)-L12) 

)-U12) 

)-L13) 

)-U13) 

(4)**2*X(9)**2)*(K11**2*K21**2  )-»-2  .  *R*G1**3 

(4)*X(5)  +  X<9)*X(m )*<K11*K12*K21*K22) 

(4)*X(6)+X(9)*X<11) )*<K11*K13*K21*K23) 

11*HH-K21*H2+X(1)*X(4)*(K11**2+K21**2)  ) 

2)*X(4)*(Kll*K12-»-K21*K22) 

3)*X(4)*(K11*K13+K21*K23) 


11*H6+K21*H7+X(1  )*X(9)*<K11**2*K21**2)  ) 
X< 2)*X<9)*(K11*K12+K21*K22) ) 
X(3)*X(9)*(K11*K13+K21*K23) ) 


(5  )**2*X(  10)**2)*(K12**2-»-K22**2)*-2.*R* 


=Y2* 
=Y2* 
=Y2* 
=Y2* 
=0.C 
=0.C 
=Y2* 
)  =  Y2 
)  =  Y2 
)=C. 
)  =  C. 
=Y2* 


<X(5)*X(6)+X< 10 )*X< 11))*(K12*K13*K22*K23) 
X<  1) *X(5)*(K11*K12+K21*K22) 
(K12*HH-K22*H2+X(2)*X(5)*XK12**2+K22**2  )  ) 
X(3)*X(5)*(K12*K13*K22*K23) 


X(1)*XU0)*(K11*K12*K21*K22) 

*(K12*H6+K2  2*H7+X(2)*X(  1^  )  *  (  «1  2**2+K22**2  )  ) 

*X<3  )*X(  1C)*(K12*K13*K22*K23) 

0 

C 

(X(6)**2+X( ll)**2)*(K13**2*K23**2)*2.*P* 


Y2*X( 1) 
Y2*X(2) 
Y2*(K13 
0.C 

c.c 

Y2*X<  1) 
=Y2*X(2 
=Y2*(K1 
=  0.C 
=C  .  C 
Y2*((  X( 
G5**3) 
Y2*X<1) 
Y2*X(  1) 
Y2*K31* 
Y2*K31* 
CO 

=c.c 
=c.c 

=  C0 

=0.0 

Y2*UX< 

G7**3) 

Y2*X(2) 

Y2*K42* 

Y2*K42* 

0.0 

=  0.0 

=c.c 

=0.0 
=  C  •  0 
Y2* ((X( 
G9**3) 
Y2*K53* 
Y2*K53* 
0.0 
=  0.0 


*X(6)*(K11*K13+K21*K23) 
*X(6)*(K12*K13+K22*K23) 
*H1  +  K23*H2«-X<3)*X(6)*(K13**2+K23**2  )  ) 


*X( 11)*( K11*K13+K21*K2  3) 
)*X(  11)*<K12*K13*K22*K23) 
3*H6+K23*H7+X(3)*X(11)*(K13**2+K 23**2) ) 


l)**2)MKll**2-t-K21**2)-t-K31**2)-»-2.*R* 


*X(2  )*<K11*K12+K21*K22) 

*X(3)*(K11*K13+K21*K23) 

K34 

K35 


2)**2)*(K12**2+K22**2)+K42**2)+2.*R* 

*X(3)*<K12*K13+K22*K23) 

K44 

K45 


3)**2)*(K13**2*K23**2)+K53**2)-»-2.*R* 

K54 
K55 
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S(6,ll)=0.0 

S(6,12)=C.C 

S(6,13  )=f  .C 

S( 7,7)=Y2*(K34**2  +  K44**2«-K54**2)*2.*R*<G10**3-G11**3) 

S(7,8)=Y2*(K34*K35+K44*K45+K54*K55) 

S(7,9)=0.C 

S( 7,10)=C.C 

s<7tii)*c.r 

S(7,12)=CC 

S<7, 13)=C.C 

S(8,8)=Y2*(K35**2*K45**2*K5  5**2)-»-2.*R1MG12**3-G13**3) 

S(8,9)=C.C 

S(8,1C  )=C.C 

S(8,ll  )=C.C 

S(8,12)=C.O 

S(8,13)=r.C 

S  (  9 ,  9  )  =  Y  2  *  ( (  X  (  1 )  *  *  2  )  *  ( K 1 1  *  *  2  +  K  2  1  *  *  2  )  *■  K  3 1  *  *  2  )  ♦  2 .  *  P  * 
$(G14**3-G15**3) 

S(9,1G)=Y2*X(1 )*X(2)*<K11*K12+K21*K22) 

S(9, 11 )=Y2*X(1 )*X(3)*(K11*K13+K21*K23) 

S( 9,12 )=Y2*K31*K34 

S(9.13  )=Y2*K31*K35 

S(  10,10) =Y2*( <  X(2)**2)*(K12**2+K22**2  H-K42  **2  I  «-2.  *R* 
£  I r 1 fe ♦ ♦ 3—  C \ 7**3 ) 

S(  lCtl  1)=Y2*X(2)*X(3)*{K12*K13+K2  2*K23) 

S(  10,12)=Y2*K42*K44 

S(  1C  ,  13)=Y2*K42*K45 

S( 11,11 )=Y2*< ( X(3)**2)*(K13**2*K23**2)+K5  3**2)«-2.*R* 
$(G18**3-G19**3) 

S(  11,12)  =  Y2*K53*K54 

S(lltl3)=Y2*K53*K55 

S<  12,12)=Y2*(K34**2  +  K44**2  +  K54**2)*2.*R*<G2r>**3-G21**3 
$) 

S( 12,13)=Y2*(K34*K35  +  K44*K45+K54*K5t>) 

S(  13,13)=Y2*(K35**2+K45**2+K5  5**2)*2.*R*(G22**3-G23**3 
$) 

DC    1     1=1, NV 

on   i  j=i,nv 

S(J,  I)=SU,J) 
1  CONTINUE 
RETURN 
END 

SUBROUTINE  MI NV( A, N,D, L ,M) 

DIMENSION  A< 1»  ,L( 1),M( 1) 

DOUBLE  PRECISION  A, 6, B IGA, HOLD 
C 

C  PURPOSE 

C  INVERT  A  MATRIX 
C 

C  USAGE 

C  CALL  MINV( A,N,D,L,M) 
C 

C  DESCRIPTION  OF  PARAMETERS 

C  A  -  INPUT  MATRIX,  DESTROYED  IN  COMPUTATION  AND 

C  REPLACED  BY  RESULTANT  INVERSE. 

C  N  -  ORDER  OF  MATRIX  A 

C  D  -  RESULTANT  DETERMINANT 

C  L  -  WORK  VECTOR  OF  LENGTH  N 

C  M  -  WORK  VECTOR  OF  LENGTH  N 
C 

C  REMARKS 

C  MATRIX  A  MUST  BE  A  GENERAL  MATRIX 
C 

C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

C  NONE 
C 

C  METHOD 

C  THE  STANDARD  GAUSS-JORDAN  METHOD  IS  USED,  THE 

C  DETERMINANT  IS  ALSO  CALCULATED.  A  DETERMINANT  OF 

C  ZERO  INDICATES  THAT  THE  MATRIX  IS  SINGULAR. 
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c 
c 

C  SEARCH    FOR    LARGEST    ELEMENT 

C 

D=1.C 

NK=-N 

DO    80    K=1,N 

NK=NK+N 

L(K)=K 

M(K)=K 

KK=NK+K 

BIGA=A(KK) 

no    20    J=K,N 

IZ=N*( J-l) 

DO    2C     I  =  K, N 

IJ=IZ+I 
10    IF(    DA6S(BIGA)-DABS(A(  UN)     15,20,2} 
15    BIGA=A(IJ) 

L(K)=I 

M(K)=J 
20    CONTINUE 
C 

C  INTERCHANGE    ROWS 

C 

J=L(K) 

IF(J-K)    35,35,25 
25    KI=K-N 

DO    30     1=1, N 

KI=KI+N 

HOLD=-A(KI  ) 

JI=KI-K+J 

A(KI  )=A(JI  ) 
30    A(  JI)     =HOLD 
C 

C  INTERCHANGE    COLUMNS 

C 

35    I=M(K) 

IF(I-K)    45,45,38 
38    JP=N*(I-1) 

DO    40    J=1,N 

JK=NK+J 

JI=JP*J 

HOLD=-A( JK) 

A( JK  )  =  A( JI  ) 
40    A(JI)    =HOLD 
C 

C  DIVIDE    COLUMN    BY    MINUS    PIVOT    (VALUE    OF    PIVOT 

C  ELEMENT     IS    CONTAINED    IN    BIGAI 

C 

45  IF(BIGA)     48,46,48 

46  0=0. C 
RETURN 

48    DO    55    1=1, N 

IF  U-K)    5C  ,55,50 
50    IK=NK+I 

A<  IK)  =  A(  IK)/(-BIGA) 
55    CONTINUE 


c 

c 
c 

REDUCE    MATRIX 

DO    65    1  =  1, N 

IK=NK+I 

HOLD=A(IK) 

I J=I-N 

DO    65    J=1,N 

IJ=IJ+N 

IF(I-K)    6C,65,60 

60 

IF(J-K)    62,65,62 

62 

KJ=IJ-I+K 

A<  IJ)=HOLD*A(KJ)*A(IJ) 

65 

CONTINUE 
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C         DIVIDE  ROW  BY  PIVOT 
C 

KJ=K-N 

DO  75  J=1,N 

KJ=KJ+N 

IF(J-K)  7C ,75, 7C 
70  A(KJ)=A(KJ J/BIGA 
75  CONTINUE 
C 

C  PRODUCT  OF  PIVOTS 

C 

D=C*BIGA 
C 

C  RFPLACE    PIVOT    BY    RECIPROCAL 

C 

A<KK)=1.C/BIGA 
80    CONTINUE 
C 
C  FINAL    ROW    AND    COLUMN     INTERCHANGE 


C 


K=N 
100    K=(K-1) 

IF(K)     15C, 150,105 
105    I=L(K) 

IF(I-K)     12C, 120,108 
108    JQ=N*( K-l) 

JR=N*< 1-1 ) 

DO    110    J=1,N 

JK=JQ+J 

HOLD=A(JK) 

JI=JR+J 

A( JK)=-A( JI) 
110    At JI  )    =HOLD 
120    J=M(K) 

IF(J-K)     ICC, 100, 125 
125    KI=K-N 

DC    13C     1=1, N 

KI=KI+N 

HOLD=A(KI) 

JI=KI-K+J 

A(KI )=-A(JII 
130    A( JI  )    =HOLD 

GO    TO    ICC 
150    RETURN 

END 
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